Marxan in R

ESM 270 Tutorial

Introduction

Designing reserve networks requires considering both ecological and socioeconomic factors. Marxan is a commonly used spatial planning tool that quantifies these factors to generate efficient reserve networks. Here, we will use species and parcel data from a 2014-15 Bren Group Project in the Morro Bay Watershed to examine the optimate regional reserve network.

Priotitizr

Prioritizr is an R package that, like Marxan, can guide systematic reserve design by solving conservation prioritization problems. In this tutorial, we will be using the marxan_problem() function which can read Marxan input data. More information on the prioritizr package can be found at https://mran.microsoft.com/snapshot/2018-03-04/web/packages/prioritizr/vignettes/prioritizr_basics.html.

A really useful tutorial in using prioritizr, specifically in using summed solutions, can be found here: https://cran.r-project.org/web/packages/prioritizr/vignettes/tasmania.html. This resource helped frame this tutorial.

Gurobi

For the solve() function, you must have a solver installed. For the the process with summing and multiple runs process, this must be gurobi. This requires following the instructions below.

1. Download gurobi (https://www.gurobi.com/downloads/gurobi-optimizer-eula/)
2. Request an academic license (https://www.gurobi.com/downloads/end-user-license-agreement-academic/)
3. Verify license in terminal by copying and pasting information in your acount > your licenses
4. Must install r package withinstall.packages('/Library/gurobi902/mac64/R/gurobi_9.0-2_R_3.6.1.tgz', repos=NULL) (if on a mac and downloaded the most recent version of gurobi, don’t download newest version of R! (4.0.0; at least as of 5/2020 this was not functional))
5. must install slam package with `install.packages(“slam”)
6. add gurobi and slam libraries

Data

We will be using a case study (and the data they prepared) from a past group project on the Morro Bay National Estuary Program.

Species data

• MorroBay_species_polyons: shapefiles of species locations
• MorroBay_species_pts: point data of species locations
• MorroBay_habitats: shapefile of location of key habitats
• morro-bay-spec.csv: information on species/habitat ID (the conservation feature), target value for that species (these values are 30% of the # of parcels containing that species or 10% for habitats), and the species name. The ‘spf’ column is the ‘species penalty factor’ and is a multiplier that determines the penalty if the conservation objective is not met
• spec-name-status.csv: information on species status

Parcel data

• morro-bay-pu.csv: includes information on the planning unit ID, the cost of that planning unit, and each parcel’s status (value 0 = available for prioritization, value 3 = locking in, 3 = locked out)
• morro-bay-puvspr.csv: includes information on the amount (presence/absence) of each species in each planning unit (pu)

Getting started

Attach packages

library(here)
library(janitor)
library(prioritizr) # marxan
library(sf) # spatial features
library(slam) # to use gurobi
library(gurobi) # solver
library(ggmap) # basemaps
library(patchwork) 
library(tidyverse)

If you do not have any of these packages installed, run the following in the console: install.packages("package name")

Note: this will not work for installing gurobi. Instructions are in the introduction section.

Read in data

All .csv files in this tutorial are unaltered from the .xlsx files provided in the R drive. They have only been renamed and saved as csvs. The “data” folder is identical to the MorroBay_data folder provided in the R drive, it has only been renamed.

Species and pu information

spec <- read_csv(here("data", "morro-bay-spec.csv")) %>% 
  head(140) %>%  # mine reads in extra blank rows - skip if yours does not 
  rename("amount" = "target") # target is called "prop" (relative) or "amount" (absolute)

pu <- read_csv(here("data", "morro-bay-pu.csv")) %>% 
  select(1:3) # mine reads in an extra blank column - skip if yours does not

puvsp <- read_csv(here("data", "morro-bay-puvspr.csv")) %>% 
  select(1:3) %>% 
  head(11849)

status <- read_csv(here("data", "spec-name-status.csv")) %>% 
  select(1:3) %>% 
  head(140)

Polygons

parcels <- read_sf(dsn = here("data"), layer = "MorroBay_parcels") %>% 
  clean_names()

ggplot(data = parcels) +
  geom_sf()

Analyze all species

Run Marxan

The function marxan_problem() in the prioritizr package gives us the “canned” approach that works for our purposes. If more customizations are desired, feel free to explore the problem() function instead.

Arguments for marxan_problem (for more information, see ?marxan_problem):

• x: planning units to use in reserve design and their cost (x = pu)
• spec: conservation features. Must contain columns for id, name, and prop/amount (spec = spec)
• puvspr: the amount of each feature in each planning unit. Must contain columns for pu, species, and amount. (puvspr = puvsp)
• bound: boundary data (we will not use this; bound = NULL)

Also note that in solver, we will use 100 runs (number_solutions = 100)

marxan_1 <- marxan_problem(x = pu,
                         spec = spec, 
                         puvspr = puvsp, 
                         bound = NULL)

marxan_1_problem <- marxan_1 %>% 
  add_gurobi_solver(gap = 0.15) %>% 
  add_pool_portfolio(method = 2, number_solutions = 100) 
# gap = optimality gap. 0.15 gives us results closest to that of using marxan software. Learn more under ?add_gurobi_solver 
# method 2 finds a specified number of solutions that are nearest to optimality. Learn more under ?add_pool_portfolio

marxan_1_soln <- solve(marxan_1_problem)

## Gurobi Optimizer version 9.0.2 build v9.0.2rc0 (mac64)
## Optimize a model with 140 rows, 670 columns and 11849 nonzeros
## Model fingerprint: 0xd2c5bccb
## Variable types: 0 continuous, 670 integer (670 binary)
## Coefficient statistics:
##   Matrix range     [1e+00, 1e+00]
##   Objective range  [2e+00, 5e+06]
##   Bounds range     [1e+00, 1e+00]
##   RHS range        [3e-01, 2e+02]
## Found heuristic solution: objective 4.238988e+07
## Presolve removed 116 rows and 205 columns
## Presolve time: 0.01s
## Presolved: 24 rows, 465 columns, 3066 nonzeros
## Variable types: 0 continuous, 465 integer (465 binary)
## Presolve removed 1 rows and 0 columns
## Presolved: 23 rows, 465 columns, 3058 nonzeros
## 
## 
## Root relaxation: objective 2.840941e+07, 14 iterations, 0.00 seconds
## 
##     Nodes    |    Current Node    |     Objective Bounds      |     Work
##  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
## 
## *    0     0               0    2.840941e+07 2.8409e+07  0.00%     -    0s
## 
## Optimal solution found at node 0 - now completing solution pool...
## 
##     Nodes    |    Current Node    |      Pool Obj. Bounds     |     Work
##              |                    |   Worst                   |
##  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
## 
##      0     0          -    0               - 2.8409e+07      -     -    0s
##      0     0          -    0               - 2.8409e+07      -     -    0s
##      0     2          -    0               - 2.8409e+07      -     -    0s
## 
## Explored 1958 nodes (383 simplex iterations) in 0.15 seconds
## Thread count was 1 (of 4 available processors)
## 
## Solution count 100: 2.84094e+07 2.84097e+07 2.842e+07 ... 3.34128e+07
## 
## Optimal solution found (tolerance 1.50e-01)
## Best objective 2.840941400000e+07, best bound 2.840941400000e+07, gap 0.0000%

Sum solutions

marxan_1_ssoln <- marxan_1_soln %>% 
  mutate(sum = rowSums(.[6:105])) %>% 
  select(id, cost, status, locked_in, locked_out, sum)

hist(marxan_1_ssoln$sum,
     main = "Selection Frequencies",
     xlab = "Number of runs that the unit was selected")

Create map

Join polygons with output

parcels_marxan_1 <- inner_join(parcels, marxan_1_ssoln, by = "id")

ggplot(data = parcels_marxan_1) +
  geom_sf(aes(fill = sum),
          color = "white",
          size = 0.05) + 
  scale_fill_gradient(low = "slategray2",
                      high = "navy") + 
  labs(fill = "Summed \nSolution") +
  theme_minimal()

Add basemap

This step is optional, but will make your map look better and is good to learn for further spatial analyses in R. Unfortunately, the options for basemaps in R are limited. Here, I have used the ggmaps package because it is the most accessible, however, it is still not very accessible compared to other R packages because you’ll need to get an API key. Instructions are here: https://cran.r-project.org/web/packages/ggmap/readme/README.html

For a ggmap cheatsheet, including the different basemaps in ggmaps, see: https://www.nceas.ucsb.edu/sites/default/files/2020-04/ggmapCheatsheet.pdf

morrobay <- get_map(location = c(lon = -120.7665, lat = 35.335),
                    zoom = 12,
                    maptype = "terrain-background", # background means no references, omit if want references
                    source = "google")

all_spec <- 
  ggmap(morrobay) +
    geom_sf(data = parcels_marxan_1,
            aes(fill = sum),
          color = "white",
          size = 0.1,
          alpha = 0.85,
          inherit.aes = FALSE) + 
  coord_sf(crs = st_crs(4326)) +
  scale_fill_gradient(low = "slategray2",
                    high = "navy") +
  labs(title = "All Species",
       fill = "Summed \nSolution",
       x = NULL,
       y = NULL) +
  theme_minimal()

all_spec

Analyze endangered species

Create dataframes

Here we create spec and puvsp data frames similar to those we had for all species, but containing only endangered or threatened species.

end_status <- status %>% 
  mutate("endangered" = case_when(
    str_detect(status, pattern = "endangered") == TRUE ~ "yes",
    str_detect(status, pattern = "threatened") == TRUE ~ "yes",
    T ~ "no")) %>% 
  filter(endangered == "yes")

end_spec <- merge(end_status, spec, by = "id") %>% 
  select(id, amount, spf, name.x) %>% 
  rename("name" = "name.x")

puvsp_id <- puvsp %>% 
  rename("id" = "species")

end_puvsp <- merge(end_status, puvsp_id, by = "id") %>% 
  rename("species" = "id")

Run Marxan

marxan_end <- marxan_problem(x = pu,
                         spec = end_spec, 
                         puvspr = end_puvsp, 
                         bound = NULL)

marxan_end_problem <- marxan_end %>% 
  add_gurobi_solver(gap = 0.15) %>% 
  add_pool_portfolio(method = 2, number_solutions = 100)

marxan_end_soln <- solve(marxan_end_problem)

## Gurobi Optimizer version 9.0.2 build v9.0.2rc0 (mac64)
## Optimize a model with 18 rows, 670 columns and 970 nonzeros
## Model fingerprint: 0x6e88f508
## Variable types: 0 continuous, 670 integer (670 binary)
## Coefficient statistics:
##   Matrix range     [1e+00, 1e+00]
##   Objective range  [2e+00, 5e+06]
##   Bounds range     [1e+00, 1e+00]
##   RHS range        [9e-01, 1e+02]
## Found heuristic solution: objective 4.174920e+07
## Presolve removed 14 rows and 203 columns
## Presolve time: 0.00s
## Presolved: 4 rows, 467 columns, 495 nonzeros
## Variable types: 0 continuous, 467 integer (467 binary)
## Found heuristic solution: objective 3.534249e+07
## Presolved: 4 rows, 467 columns, 495 nonzeros
## 
## 
## Root relaxation: objective 2.677119e+07, 1 iterations, 0.00 seconds
## 
##     Nodes    |    Current Node    |     Objective Bounds      |     Work
##  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
## 
## *    0     0               0    2.677119e+07 2.6771e+07  0.00%     -    0s
## 
## Optimal solution found at node 0 - now completing solution pool...
## 
##     Nodes    |    Current Node    |      Pool Obj. Bounds     |     Work
##              |                    |   Worst                   |
##  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
## 
##      0     0          -    0               - 2.6771e+07      -     -    0s
##      0     0          -    0               - 2.6771e+07      -     -    0s
##      0     2          -    0               - 2.6771e+07      -     -    0s
## 
## Explored 4914 nodes (1716 simplex iterations) in 0.41 seconds
## Thread count was 1 (of 4 available processors)
## 
## Solution count 100: 2.67712e+07 2.67732e+07 2.67736e+07 ... 3.14836e+07
## 
## Optimal solution found (tolerance 1.50e-01)
## Best objective 2.677119300000e+07, best bound 2.677119300000e+07, gap 0.0000%

Sum solutions

marxan_end_ssoln <- marxan_end_soln %>% 
  mutate(sum = rowSums(.[6:105])) %>% 
  select(id, cost, status, locked_in, locked_out, sum)

hist(marxan_end_ssoln$sum,
     main = "Selection frequencies",
     xlab = "Number of runs that the units were selected")

Create map

Join polygons with output

parcels_marxan_end <- inner_join(parcels, marxan_end_ssoln, by = "id")

ggplot(data = parcels_marxan_end) +
  geom_sf(aes(fill = sum),
          color = "white",
          size = 0.05) + 
  scale_fill_gradient(low = "slategray2",
                      high = "navy") + 
  labs(fill = "Summed \nSolution") +
  theme_minimal()

Add basemap

end_spec <- 
  ggmap(morrobay) +
    geom_sf(data = parcels_marxan_end,
            aes(fill = sum),
          color = "white",
          size = 0.1,
          alpha = 0.85,
          inherit.aes = FALSE) + 
  coord_sf(crs = st_crs(4326)) +
  scale_fill_gradient(low = "slategray2",
                      high = "navy") + 
  labs(title = "Endangered & Threatened Species",
       fill = "Summed \nSolution",
       x = NULL,
       y = NULL) +
  theme_minimal()

end_spec

Combine maps with patchwork

This step isn’t necessary, but is really convenient if you want to present your maps together in your report. patchwork uses PEMDAS to combine ggplots together into one graphic. For example, plot_1 + plot_2 will make an image of the plots side by side, whereas plot_1 / plot_2 will make an image of plot_1 over plot_2. For more information, see https://github.com/thomasp85/patchwork.

First, make the all_spec graph without a legend so that the combined image has only one legend.

all_no_lgnd <- all_spec +
  theme(legend.position = "none")

Next, combine the graphs with patchwork

spec_graphs <- all_no_lgnd + end_spec

spec_graphs

Finally, save the image to add to your report!

ggsave("spec-graphs.png")